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Abstract. We study quantum entanglement in the ground state of the Affleck- 
Kennedy-Lieb-Tasaki (AKLT) model defined on two-dimensional graphs with reflection 
and/or inversion symmetry. The ground state of this spin model is known as the 
valence-bond-solid state. We investigate the properties of reduced density matrix of a 
subsystem which is a mirror image of the other one. Thanks to the reflection symmetry, 
the eigenvalues of the reduced density matrix can be obtained by numerically 
diagonalizing a real symmetric matrix whose elements are calculated by Monte Carlo 
integration. We calculate the von Neumann entropy of the reduced density matrix. The 
obtained results indicate that there is some deviation from the naive expectation that 
the von Neumann entropy per valence bond on the boundary between the subsystems is 
In 2. This deviation is interpreted in terms of the hidden spin chain along the boundary 
between the subsystems. In some cases where graphs are on ladders, the numerical 
results are analytically or algebraically conflrmed. 
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1. Introduction 

In recent years, the study of entanglement in quantum many-body systems lias become 
a common issue in condensed matter physics, statistical mechanics, and quantum 
information theory. Entanglement is a purely quantum mechanical phenomenon in 
which quantum systems are linked together even when they are spatially separated so 
that one system cannot be correctly described without full mention of its counterpart. 
This concept was first introduced by Schrodinger in [1]. From the viewpoint of quantum 
information science, entanglement is a fundamental measure of how much quantum 
effects we can observe and use to control one quantum system by another, and it is 
the primary resource for quantum information processing and communication [21 |3l H] . 
On the other hand, in condensed matter physics and statistical mechanics, the 
concept of quantum entanglement has recently been used to investigate quantum phase 
transitions [5] (see [6] for a review), topological order [71|8l|9], and macroscopic properties 
of solids [in]. A fundamental and practical problem common in all the areas is how to 
detect entanglement [11] and quantify the degree of it. A general approach is to study 
the reduced density matrix of a certain subsystem in an entangled state. The concept 
of the reduced density matrix was first introduced by Dirac in [12]. There are many 
kinds of measures of entanglement related to the reduced density matrix (see [13] and 
references therein). Among them, one of the most popular characteristic functions is 
the von Neumann entropy (entanglement entropy) of a reduced density matrix [HI [15] . 
It serves as the measure of entanglement for a pure state. 

Generally speaking, to elucidate quantum entanglement in many-body systems is 
a formidable task since we need a precise description of many-body quantum states. So 
far, much insight has been obtained by studying exactly solvable models in which there 
is a possibility to obtain the reduced density matrix and calculate the von Neumann 
entropy exactly or approximately. Representative models are the harmonic models 
[Ml [IT], the XY spin chain [HI [13 [201 [2ll [221 [23], the XXZ spin chain [21 [23 [26], 
the Calogero- Sutherland model [27], one- dimensional critical models [28], conformal 
field theories [291 [30], the Affleck-Kennedy-Lieb-Tasaki (AKLT) model [SH [321 [33] , the 
Kitaev model [M], and the quantum dimer (Rokhsar-Kivelson) model [35]. Among 
them, the AKLT model [361 EZ] has some important meanings since this model is 
directly related to one of the schemes of quantum computation, measurement based 
quantum computation [381 EH] . Historically, the AKLT model and its exact ground state 
known as the valence-bond- solid (VBS) state were proposed to understand ground-state 
properties of the Haldane gap systems [40l[4l]. The authors (AKLT) rigorously proved 
that the VBS state is a unique ground state and there is an energy gap immediately 
above the ground state. Exponential decay of correlation functions in the VBS ground 
state has also been shown. There is an intimate relation between the VBS states and 
fractional quantum Hall states which has been revealed by using a spin coherent state 
representation [42]. A hidden topological order in the VBS state is characterized by 
string order parameters [l3] and existence of edge states which are degenerate ground 
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states in a chain with open boundary conditions US]. The VBS picture and the 
edge states were also experimentally supported by electron spin resonance (ESR) in the 
S = 1 Haldane chain Ni(C2H8N2)2N02(C104) (NENP) 

Although most of the exactly solvable models are limited to one dimension, 
the AKLT model and VBS ground state can be generalized to higher dimensions 
and/or inhomogeneous (non-translational invariant) lattices [371 SHI SHI [50], [51]. This 
feature is especially important in connection with the measurement based quantum 
computation as we discussed later. For studying nature of the VBS states, the non- 
trivial observation has been done by Kliimper et al, through the g-deformation of the 
AKLT model [52], |53l [M] . They found that the one-dimensional (Id) spin-1 VBS state 
can be expressed in a form of the matrix product state (MPS). This representation 
gives an efficient method to calculate correlation functions. Mathematical theory and 
generalization of the VBS states, i.e., finitely correlated states (FCS), were essentially 
developed by Pannes et al, in [55] [56] [57] [58] . Another interesting generalization of the 
VBS state is to replace SU(2) spins with other Lie algebra generators. In this direction, 
there have been constructed many versions of the AKLT model such as SU(3) [59] 
and SU(n) modes [HDl [EH [62], SU(2n)-extended model [63], symplectic model [6^ . 
S0(5) model [65], S0(2S' -|- 1) model [661 [67], supersymmetric model [68], and several 
g-deformations [69] [701 [HI [72] [73] [71] . 

The VBS state has recently been attracting renewed interest from the viewpoint 
of quantum information theory. In particular, it has been proposed that the VBS state 
can be used as a resource state in measurement-based quantum computation (MQC). 
MQC is one of the schemes of quantum computation in which we can perform universal 
quantum computation using only measurements as computational steps (see [75] for 
an introduction). There are two kinds of measurement based models. One is the 
teleportation quantum computation (TQC) [TB] [77] [78]. The other is the so-called 
one-way quantum computation (IWQC) or cluster state computation [TQ] [HDl El 182] . 
In this model, any quantum gate array can be implemented as a pattern of single- 
qubit measurement on a highly entangled cluster state. Note that IWQC with one- 
dimensional cluster states is insufficient for universal quantum computation and thus 
two-dimensional cluster states are inevitably needed. It has then been proved in [3H] 
by Verstraete and Cirac that TQC and IWQC are equivalent using the VBS picture. 
Gross, Eisert and collaborators have introduced novel schemes for MQC based on finitely 
correlated or tensor network state beyond the cluster state [83] [HI]. Another recent 
important development by Brennen and Miyake is that the VBS state (more precisely 
dynamically coupled AKLT chains) can be used as a resource state in MQC instead 
of the cluster state [39]. The VBS state is more advantageous to IWQC than the 
cluster state since it is robust against noise due to the energy gap. Implementations of 
the AKLT Hamiltonian using bosonic atoms [HS] [HS] or spin-1 polar molecules [H7] in 
optical lattices have also been proposed. 

In this paper, motivated by considerable current interest from several viewpoints, 
we study entanglement properties of VBS states on two-dimensional (2d) graphs which 
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is relevant to MQC Before discussing the 2d VBS states, let us briefly summarize 
the obtained results in Id VBS states. The reduced density matrix and entanglement 
properties in the Id spin-1 VBS state was flrst studied in [3l]. The authors found that 
the von Neumann entropy of the block with the rest of the chain approaches a constant 
value exponentially fast, as the size of the block increases. The constant value is given 
by 2 In 2. Then, the Id integer-S* VBS states were studied in [33]. In this case, the 
von Neumann entropy of the large block with the rest is given by 2 ln(5' + 1) and this 
value is interpreted in terms of the edge states, the degenerate ground states of the 
open AKLT chain. This picture leads to an idea that the reduced density matrix of 
the block in the AKLT chain is exactly spanned by the degenerate ground states of the 
block Hamiltonian (see Eq. ( 127|) for the deflnition). This idea has been conflrmed in the 
series of papers [62l [HHl [HH] and was summarized in the review [90]. Another interesting 
quantities such as the entanglement length [91], the geometric entanglement [92], the 
single-copy entanglement [93], and the boundary effect [91] in the Id spin-1 VBS state 
have also been studied. More generally, the entanglement of formation in the Id FCS 
was estimated in [95] and an area law in Id gapped quantum system was proven in |96j . 
Compared to the Id VBS states, few results have been obtained in 2d VBS states so far. 
Entanglement properties of the VBS state on the Cayley tree [371 [97] was studied in [98] . 
The authors showed that the von Neumann entropy of the block does not depend on 
the whole size of the system and its asymptotic value is linearly proportional to the 
number of valence bonds crossing the boundary. This gives an explicit proof of the area 
law in this system (see [99] for a review of area laws and the entanglement entropy). 
General entanglement properties of the VBS state on an arbitrary graph [371 IHl EQ] have 
been studied in [lOOj . The authors showed that the eigenspace of the reduced density 
matrix of the block is spanned by the degenerate ground states of the block Hamiltonian. 
The authors have shown that the von Neumann entropy of a large block with the rest 
approaches a value less than \n{deg.), where (deg.) is the ground-state degeneracy of 
the block Hamiltonian. Although some general results for 2d VBS states have been 
obtained so far, a quantitative analysis of entanglement was missing since a powerful 
analytical method such as a transfer matrix technique does not work sufficiently in the 
analysis of 2d systems. Therefore, entanglement properties of the 2d VBS states must 
be analyzed with the aid of numerics. 

In this paper, we study entanglement properties of the 2d VBS states on symmetric 
graphs. Symmetries such as reflection and inversion enable us to develop an efficient 
method to study the reduced density matrix. Due to the symmetry, one can find the 
relation between the eigenvalues of the reduced density matrix and those of the overlap 
matrix (see Eqs. ( !T5llT6|) for the deflnition). The overlap matrix is a real symmetric 
and its elements are given by the inner product of the degenerate ground states of the 
block Hamiltonian which are linearly independent but may not be orthonormal. The 
inner product can be, for example, calculated by Monte Carlo integration using the 
Schwinger boson representation of the VBS state. Even without numerical integration, 
one can formally write the overlap matrix as a matrix-valued correlation function in 
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the VBS state restricted to one of the subsystems. From this, we conjecture that the 
von Neumann entropy per valence bond (boundary site) is strictly less than In 2 even 
in the 2d infinite system. A holographic interpretation of the deviation from In 2 is 
given in terms of the hidden spin chain along the boundary between two subsystems 
(refiection axis). To confirm this conjecture, we numerically study the eigenvalues 
of the overlap matrix by exact diagonalization. Since the dimension of the overlap 
matrix does not depend on the size of the whole system but is given by (deg.), we can 
study relatively large systems. From the eigenvalues of the reduced density matrix, we 
calculate the von Neumann entropy. We find that the von Neumann entropy per valence 
bond on the boundary is strictly less than In 2 for finite size systems. We also find a 
function which fits the obtained numerical results well. The long-distance behavior of 
this function strongly supports our conjecture. These numerical results are analytically 
or algebraically confirmed for the graphs with a ladder geometry, where transfer matrix 
technique can be applied. 

The organization of this work is as follows. In Section 2, we review the construction 
of the AKLT model on an arbitrary graph and introduce the Schwinger boson 
representation of the VBS state. In section 3, we provide a general argument on the 
Schmidt decomposition and discuss a relation between the reduced density matrix and 
the overlap matrix. We show how to apply this method to study the von Neumann 
entropy of the VBS state. Then we make a conjecture on the von Neumann entropy in 
the 2d infinite system. In Section 4, numerical results of the VBS states on the square 
and hexagonal lattices are shown. In Section 5, the obtained numerical results are partly 
confirmed analytically for graphs in a ladder geometry. Conclusions are given in the 
last section. In Appendix, we provide a graphical interpretation of the overlap matrix 
in terms of closed loops and open strands. 

2. The basic AKLT model 

Let us first define the basic AKLT model on a connected graph. A graph consists of 
two types of elements, namely vertices (sites) and edges (bonds). As shown in Fig. [H 
a vertex is drawn as a (large) circle and an edge is as a solid line. A graph is called 
connected if every pair of distinct vertices in the graph can be connected through some 
path. Let G = {V, E) be a connected graph, where V and E are the sets of vertices and 
edges, respectively. Henceforth, we assume \V\ > 1, otherwise there is no interaction 
between spins. The spin operator Sk is located at vertex k and its spin value is denoted 
by Sk- In the basic model, we require that Sk = \zki where Zk is the coordination 
number, i.e., the number of incident edges connected to k. To guarantee uniqueness 
of the ground state, this relation must be true for any vertex /c, including boundaries. 
Let us now define the spin Hamiltonian for the basic AKLT model. The Hamiltonian is 
written as a sum of interactions on all the edges: 




(1) 



{k,l)eE 
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Here, A{k,l) is an arbitrary positive real coefficient (it may depend on the edge {k,l)) 
and the operator 7is^+Si{k, I) projects the total spin on the edge {k, I), Jk,i = Sk + Si, on 
the subspace with the highest possible spin value J^a = Sk + 5*/. Since the Hamiltonian 
in Eq. ([T]) is a linear combination of projection operators with positive coefficients, H 
is positive semi-definite. The explicit form of the projector ns^+SiikJ) in terms of Sk 
and Si is given by 

^ (kJ^- n {Sk + SiY-j{j + i) 

^s,^s,{kJ)- 11 iSk + Si)iSk + Si + l)-jij + iy 

\Sk-S,\<j<S^:+Sl-l " ^'^ ) J\J I 

Note that one can expand {Sk + SiY = 2Sk ■ Si + Sk{Sk + 1) + Si{Si + 1) and hence the 
projector defined in Eq. (|2]) is a polynomial in the scalar product {Sk ■ Si) of the degree 
S^min, where 5*111111 = Sk is the minimum of Sk and Si. As an example, consider a case 
where Sk = Si = 1, which corresponds to the S* = 1 homogeneous AKLT chain. The 
projector 7i2{k,l) is written as 

MkJ) = l{Sk-Siy + ^{Sk-Si) + ^. (3) 

The projectors Ti^{k,l) and TTi{k,l) are related to the Hamiltonians on the hexagonal 

and the square lattices, respectively: 

1 ^ ^ 2Q ^ ^ 27 ^ ^ 11 
-3(*, = ^(S. . S,)» + -(S. . S,f + _(S. . S,) + -. (4) 

In general, explicit expressions in terms of (5*^ • Si) become more complicated. 
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Let us now discuss the ground state of H in Eq. ([T]) with condition 

Sk = ^Zk. (6) 

The ground state is unique and known as the valence-bond-sohd (VBS) state. Although 
there are several possible representations for the VBS state (see [901 HOOj )- the most 
convenient one for our purposes is the Schwinger boson representation [l2l HHj |50] . In 
this representation, we introduce a pair of bosonic creation and annihilation operators 
at each vertex to realize SU(2) Lie algebra. We define a pair of bosonic operators ak 
and bk for each vertex k with the canonical commutation relations: 

[afc,af] = [bk,bl] = Sm (7) 
with all other commutators vanishing: 

[ttk, ai] = [bk, bi] = [ak, k] = [ak, b]] = 0, ^(/c, /). (8) 
Spin operators are represented by the Schwinger bosons as 

St = <^lbk, S^ = biak, S'k = li4<^k-blbk). (9) 

To reproduce the dimension of the spin-S'^ Hilbert space at each vertex k, the following 
constraint on the total boson occupation number is required: 

4«fc + bibk = 2Sk. (10) 

With this constraint in mind, the VBS ground state is written as 

|VBS)= n (4&!-&W)|vac), (11) 

where the product runs over all edges and the vacuum |vac) is defined as the direct 
product of vacuums of each vertex, i.e., 

|vac) = (S?)|vac['=l). (12) 



Note that jvact'^l) is destroyed by any annihilation operators ak or bk- Using the 
Schwinger boson representation, one can generalize the VBS state which is the ground 
state of the generalized AKLT model. By associating a positive integer Mk^i to each 
edge {k, I) of G, the generalized VBS state is written as 

|VBSGen)= n (13) 

{k,l)eE 

For the construction of the generalized AKLT Hamiltonian and the condition of the 
uniqueness of the ground state, please see original references [HI [50|, [901 HOO] . In this 
paper, we shall focus on the basic model, i.e., Mk^i = 1 for any {k, I). 

Let us now prove that the VBS state in Eq. (fTT!) is the zero-energy groud state of 
the basic AKLT Hamiltonian H. To prove this, we have only to show for any vertex k 
and edge {k, I) (i) the total power of a^ and 6^, is 25*^ so that the spin value at the vertex 
k Sk, and (ii) the maximum value of the total spin of the edge {k, I) is Sk + Si — 1. (i) can 
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be shown by expanding r.h.s. of Eq. fllip and finding tliat tfie total power of and 6^ is 
Zk = 25^. Let us prove (ii). By counting the power of a^'s minus 6^'s on the bond {k, I), 
we find that the maximal eigenvalue of + S'f is S*?; + S'; — 1, i.e., J^i < Sk + Si — 1. 
Since the state |VBS) is invariant under global rotations, the maximum value of the 
total spin of the edge {k,l) is Sk + Si — 1. Therefore, from (i) and (ii), the VBS state 
|VBS) in Eq. ffTTl) is a zero-energy ground state of the basic Hamiltonian H. It was 
shown in [HI HI] that this ground state is unique. 



3. Schmidt decomposition and VBS states on symmetric graphs 

In this section, we apply the Schmidt decomposition method to the VBS state on a 
reflection symmetric graph. We will see that an upper bound on the von Neumann 
entropy can be easily obtained by this method. We shall start with the Schmidt 
decomposition for a general state. We follow the approach by Shi et al. in [lOlj . 
Suppose that the state of a total system is written as 



a 

where A and B denote subsystems A and -B, respectively. We call the above 
decomposition pre-Schmidt decomposition. Note that the sets of states and 
linearly independent but may not be orthonormal. In this sense, the above 
expression is not of the form of the Schmidt decomposition. To Schmidt decompose the 
above state, we now define overlap matrices M^"^! and M^^l as 

{M^^^U^{4>^^^\<pf), (15) 

(M[^'W = (0L'^^l0!f'), (16) 
respectively. Their spectral decompositions are given by 

Ml^l = VDl^lvt, {X^X = XX^ = 1), (17) 

Ml^l = YD^^^Y\ {Y^Y = YY^ = 1), (18) 

where D'^^'^ and D^^^ are diagonal matrices, i.e., {D^"'^)rT' = Srr'd^r^ {a = A or B). Using 
X and V, one can obtain the orthonormal sets in A and B as 

Then, Eq. f ll4p can be written as 



I*) = E ^d^rdf\x^YU\er) ® I/,), (20) 

TTI 

where denotes matrix transpose. So far, we have considered a general case. Let us now 
focus on the special situation where M^^^ and M^^"^ are the same, i.e., 

^[A] ^ j^lB] ^ 

Furthermore, if M is a real symmetric matrix, it can be diagonalized by some orthogonal 
matrix O as M = ODO^ . In such Eq. fl20p becomes a simpler form: 

\^)=Y,dr\er)®\fr). (21) 
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Figure 2. VBS state on a reflection symmtric graph. and As denote the sets of 
sites (vertices) on the boundaries of subgraphs A and B. Any site i € A^ has a unique 
reflection symmetric partner i & B. 



Note that = rfi^' = SfK One may think that the above situation is not generic. 
However, as we will see later, it is enough for our purpose to restrict ourselves to such a 
special case. Let us now consider the density matrix for the total system consisting of 
A and B. It is defined by p = |\E')(\E'|/(\E'|\E'). Then the reduced density matrix for A is 
given by = tr^ p. In our case, it is given by 

Pa = '^{fvlplfv) 
_ Er dl\er){er\ 

Therefore, the von Neumann entropy in this bipartite partitioning is given by 

5 = -trp^lnpA = -^Prlnpr, with = (23) 

Therefore, we can obtain the von Neumann entropy in terms of the overlap matrix M. 
To obtain 5, we need all the eigenvalues of M. 

Let us now apply the Schmidt decomposition method to study entanglement in VBS 
states on graphs with symmetries such as reflection and/or inversion. For simplicity, we 
shall focus on the graphs with reflection symmetry. However, one can easily generalize 
our argument to graphs with inversion symmetry. We define reflection symmetry such 
that any site (vertex) in a subgraph A has a reflection partner in B and vice versa. A 
graphical example is shown in Fig. [21 In the Schwinger boson language, the VBS state 
on a reflection symmetric graph can be written as 

ivBS)= n K^-^H) n nK^-^H)i^ac)' (24) 
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where Ba and Bb are sets of edges (bonds) in subsystems A and B while Aa and 
are sets of vertices (sites) on the boundaries of A and B. We note that we only consider 
the basic VBS state which is the ground state of the basic AKLT model. In this state, 
there is one valence bond on any edge The generalization of our argument to 

generic VBS states is straightforward. To apply the Schmidt decomposition, we slightly 
change the above expression by a local gauge transformation: Uj — )■ —bj and bj — )■ aj 
for ^jEB. Then Eq. flMl) becomes 

ivBS)= n n K^-^H) n(«H+^H)i^^^)- (25) 

Since the graph is reflection symmetric, i G A^i necessarily has a reflection symmetric 
partner i & Kb (see Fig. |2]). For simplicity, we henceforth assume that any i is nearest 
neighbor to i in the graph G. Therefore, the above expression can be written as follows: 

{a} 

where {a} = {ai, 02, a|A^|} («» = ±1/2) corresponds to the spin state of the 
boundary state, i.e., degenerate ground state of the block Hamiltonian defined within 
the subsystem A 01 B. The block Hamiltonians are explicitly defined as 

Ha= Mk,l)7is,+sXkJ), Hb= MkJ)7is,+s,ik,l), (27) 

{k,l)eBA {k,l)eBB 

where 7rSk+Si{k,l) is defined in Eq. ([2]). The states and are given by 

= n n K^-^H)ivac[^^ 

= n(4)'/^^"H&I)^/^-"^ n (al&l-&H)|vac[^l), (28) 

respectively, where Ivac'"^!) dvac'^')) is the vacuum for bosons in A (B) and |vac) = 
Ivacl'^l) ® Ivac^'^l). Equations in (1251) immediately yield that overlap matrices M^^l and 
are the same. The matrix element of M(= M^^l = M'-^l) is given by 

M|,|,l^} = (vac[^]| n n(«^)'^''""H&fc)'/'""^ 

{i,j)&BA keAA 

X n (4)^/^^^H&l)^/^-^^ n («H-^M)|vac[^l)- (29) 

k&^A {i,j)&BA 

The rank of M is 21^^', where we denote the number of elements in a set S by IS*]. 
The matrix element of M is real, i.e., {M^a},{i3})* = ^{o},{/3}; because the commutations 
of bosons do not produce any complex phases. Therefore, one can easily show the 
symmetric property M{a},{i3} = M{^},{^}. 
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We now rewrite the above expression in terms of the spherical angles. This can be 
done by introducing spin coherent state representation. Let us first introduce spinor 
coordinates {uk,Vk) at each vertex (site) k: 

«, = e^<^^/2cos|, v, = e-^^^/'siJ-^, (30) 

where < 6 < n and < </> < 27r. Then for a point = 
(sin 9k cos (j)k, sin 9k sin 0^, cos 9k) on the unit sphere, the spin-S'^ coherent state is defined 
as follows: 

l^k) ^ ^^^^^i±M)!!^-|vacW). (31) 

Here we have fixed the overall phase (U(l) gauge degree of freedom) since it has no 
physical content. The coherent state is not orthogonal but complete and the resolution 
of identity is given by 

hsk+i= XI \Sk,'m){Sk,rn\ = '^^^^^ I dVLk\VLk){VLk\, (32) 

where |5'fc,m) denotes the simultaneous eigenstate of 5*1 and S^, and 125^+1 is the 
(25*^ + l)-dimensional identity matrix. By inserting this resolution of the identity, Eq. 
(!29|) can be recast in an integral form: 




X 



n . (33) 

where we have used the fact {Ya.c\a^''~%^'''^^\VLk) = \/ {^Sk)\u^k~^v^^^ . Therefore, if we 
can obtain the matrix element of M by some method such as Monte Carlo integration, 
the only thing to do is to diagonalize the matrix M. Once the eigenvalues of M, i.e., dr 
are obtained, the von Neumann entropy can be calculated through 
Eq. ( 1231) . Let us now estimate the upper bound on the von Neumann entropy. Suppose 
that all the eigenvalues of M are the same. This means that the subsystems A and B 
are maximally entangled. Then, the von Neumann entropy is given by |Aa| In 2. This 
value gives an upper bound on the von Neumann entropy, i.e., S < |A^| In 2. However, 
as we will see in the next section, numerical results indicate a deviation from this naive 
estimate and the von Neumann entropy S is strictly less than \Aa\ In 2 in 2d VBS. 

It is interesting to note that there is a holographic interpretation of the overlap 
matrix M. The matrix M can be regarded as an operator acting on the auxiliary 
one dimensional chain which is along the boundary of A, i.e., A^. From the fact that 
Clk = (sin 9k cos (pk, sin 9k sin (pk, cos 9k), one can rewrite M as 



M= n^i^« n n . (34) 
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Here we have used the fact that M'^ = M. Note that the Pauh matrices (a = x, y, z) 
act on the auxihary space, {| t'*')' I i'*^)}- Then using the following relation: 



(5', + l)(2^, + l) 



dQiQi\Qi){Qi\ = Si, {a = x,y,z), (35) 



the overlap matrix can be expressed as a matrix-valued correlation function: 

M = (VBS[^]| l[(s, + l + S,- a) iVBSt^l), (36) 

where 

iVBSt^l)^ n -&!«]) Ivad^^) (37) 

Note that the matrix rank of Si for i G is 25*4 — 1 while that for z G A \ A a is 
25*4 + 1. One can interpret the overlap matrix as a Hamiltonian acting on cr-spin space. 
The matrix element of this Hamiltonian is governed by the boundary correlation in the 
half-VBS state given in Eq. ( 13 7p . The naive upper bound value of the von Neumann 
entropy, |Ayi|ln2, corresponds to the strongly disordered limit where any boundary 
correlation among S'-spins vanishes. However, as shown numerically in |102[ 1103] , the 
spin-spin correlation functions at small distances are nonzero in the VBS-type state 
even though they decay exponentially fast. Therefore, the leading term in Eq. (136!) 
is given by a constant plus a term proportional to the spin-| Heisenberg Hamiltonian 
whose coupling J is given by the nearest neighbor S-spin correlation function. The 
eigenvalues of this matrix can be different and hence the entanglement spectrum may 
not be flat. Therefore, it is plausible that the von Neumann entropy is strictly less than 
|A^| In 2 even in the 2d infinite system. To support this conjecture, we perform numerical 
and analytical calculations for the VBS states on square and hexagonal lattices in the 
following sections. 



4. Numerical analysis of 2d VBS states 

In the previous section, we derived the integral formula for the overlap matrix (Eq. 
( 133|) ) which helps us to reduce computational cost. Let us explain this point in more 
detail. Suppose that the graph G has \V\ vertices and every vertex has the same 
spin S. Then, the needed dimension for representing the VBS state on G is (25* + 1)'^' 
without consideration of symmetry. However, if we focus on the von Neumann entropy of 
reflection symmetric VBS state, the needed dimension is greatly reduced from (25 + 1)'^' 
to 2l^'*l, where |A^| denotes the number of sites on the boundary between A and B. 
Therefore, we can study the von Neumann entropy of 2d VBS state on a relatively 
large system by combining Monte Carlo integration and exact diagonalization. After 
obtaining the overlap matrix and all its eigenvalues, the von Neumann entropy can be 
calculated according to Eq. ( l23l) . In this section, we calculate von Neumann entropies 
of VBS states on square and hexagonal lattices with open boundary conditions using 
Monte Carlo integration. 
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Figure 3. Square lattice with open boundary conditions. A^x and Ny denote the 
numbers of sites along x- and y-axes in A, respectively. 



4-.1. von Neumann entropy of VBS state on square lattice 

We first study the von Neumann entropy of the VBS state on square lattices with open 
boundary conditions. The reflection axis of the square lattice is taken to be the boundary 
between A and B as shown in Fig. [3l in Fig. |3] is the number of sites along x-axis 
while Ny is that along y-axis. For square lattice, the number of boundary sites |A^| is 
the same as Ny, i.e., |Aa| = Ny. We first study the A^a;-dependence of the von Neumann 
entropy. Figure S] shows the von Neumann entropy per boundary site, i.e., valence bond, 
for various Ny as a function of A^^;. The von Neumann entropy per valence bond, iS/|Ayi|, 
rapidly converges to certain values depending on Ny. The case of A^y = 1 reduces to 
the linear AKLT chain along x-axis. In this case, iS/|A^| = S monotonically decreases 
with increasing A^; and approaches to In 2 (= 0.6931472...) exponentially fast. This 
correctly reproduces the exact result in [31]. The obtained results show that 5/|A^| for 
Ny = 2 and Ny = 3 approach to 0.6494348 and 0.6315983, respectively. These values 
are consistent with the analytical results shown in the next section. Next we consider 
the A^y- dependence of the von Neumann entropy. Figure |5] shows the von Neumann 
entropy per boundary site as a function of Ny. We find that iS/|A^| is well fitted by the 
following function: 

where C, A and a are fitting parameters. The constant term a means the von Neumann 
entropy per valence bond in the limit of |A^| — > oo. Note that C, A and a depend on 
Nx. The blue curves in Fig. [5] represent the fitting curves. The obtained numerical 
data show good agreement with the power law behavior assumed in Eq. f l38|) . Table 
[1] shows the fitting parameters C, A and a for A^; = 1,2,3,4 and 5. As A^^; increases, 
the coefficient C increases and A decreases. The non-universal coefficient for the area 
law, a, decreases very slowly when N^ increases. Suppose that A does not vanish in the 
large Nr^ limit. Then, we can rewrite Eq. fl38|) as 

S = a\AA\+C\AA\^~^. (39) 
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Figure 4. The von Neumann entropy per valence bond as a function of for the 
square VBS states with various Ny. The dotted hues indicate the exact resuhs in 
Nx — >■ oo hmit. 
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Figure 5. Entanglement entropy per valence bond as a function of A'y for the square 
lattice VBS states with various N^. Blue curves represent the fitting curves (Eq. (|38|) ). 



The first term simply means the area law while the second term is a tiny correction 
from it. Our numerical results indicate the coefficient a is strictly less than In 2 in the 
Nx — )■ oo limit. 
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C 


A 


a 




0.0821(4) 


0.90(1) 


0.6110(4) 


iV, = 2 


0.104(1) 


0.78(2) 


0.589(1) 


iVx = 3 


0.108(1) 


0.75(2) 


0.584(1) 


iV, =4 


0.110(2) 


0.73(3) 


0.582(2) 


iV, = 5 


0.112(1) 


0.71(2) 


0.580(1) 



Table 1. The coefficients of the fitting function Eq. 



for square lattices. 



B 



Figure 6. Hexagonal lattice with open boundary condition. Nx and Ny denote the 
numbers of the sites along x- and y- axes in A, respectively. 



4-2. von Neumann entropy of VBS state on hexagonal lattice 

In this subsection, we consider the von Neumann entropy of the VBS state on hexagonal 
lattices with open boundary conditions. Similarly to the case of the square lattice, 
subsystems A and B are partitioned by the reflection axis shown in Fig. |6l in Fig. 
|6] is the number of sites along x-axis while A'j^ is that along y-axis. Since we focus on 
even-Ny cases in the present study, the number of the boundary sites is |Aa| = Ny/2. 
We first study the A^^^- dependence of von Neumann entropy. The results of the von 
Neumann entropy per boundary site, i.e., valence bond, are shown in Fig. [71 The case 
of \^a\ = 1 reduces to the S = 1 AKLT chain and reproduces the exact result in [3T| . 
The von Neumann entropy per valence bond, i5/|Aa|, for Ny = 2 and 3 approach to 
0.6890926 and 0.6876522, respectively. These values also consistent with the analytical 
results shown in Section |5l Next we study the A^y- dependence of the von Neumann 
entropy. The obtained results are shown in Fig. |H1 Those results are again well fitted 
by the scaling function given by Eq. ( 138|) . Those fitting curves are shown by the blue 
curves in Fig. [HI The fitting parameters C, A, and a for N^ = 1,2,3,4, and 5 are 
summarized in Table [2j As N^ increases, C increases and A and a decreases. Assuming 
A is nonvanishing in the large N^ limit, iS/|A^| is strictly less than In 2 even in the 
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Figure 7. Entanglement entropy per valence bond as a function of for the 
hexagonal VBS states with various |A^|. The dotted lines indicate the exact results in 
Nx — )■ oo limit. 
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Figure 8. The von Neumann entropy per valence bond as a function of Ny for the 
hexagonal VBS states with various N^- Blue curves represent the fitting curves. 



infinite 2d system. Therefore, our numerical results again supports the conjecture given 
in Sec. [3l 

5. Algebraic analysis of VBS ladders 

In this section, we shall consider the VBS states on various ladders and study the von 
Neumann entropy via another way. It is convenient to introduce the following function 
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C 


A 


a 


iV, = 1 


0.00812(1) 


0.974(4) 


0.68502(1) 


N. = 2 


0.00849(4) 


0.94(1) 


0.68465(5) 


iVx = 3 


0.00870(3) 


0.904(7) 


0.68443(3) 


iV, = 4 


0.0092(2) 


0.82(3) 


0.6838(2) 


iVx = 5 


0.00941(7) 


0.81(1) 


0.68373(7) 



Table 2. The coefficients of the fitting function Eq. 



for hexagonal lattices. 
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Figure 9. VBS state on a 2-lcg ladder. It is cut by the reflection line indicated by 
the broken line. 



of {Si}: 



'|Aa|. 




An 



Ylil + dk-s,) n (l-^i-^j), (40) 

where Sj's are classical or operator- valued vectors depending on the context. The overlap 
matrix given in Eq. (!34|) is proportional to Z{ai, ...,a\\^\). Since the only difference 
between Z and M is the overall factor, the eigenvalues of the reduced density matrix 
equal to those of the following matrix: 
Z(ai,...,a|A^l) 

Pa - —^7^ 75 (41) 

trZ(ai, ...,cr|A^|) 

where the trace is taken over the cx-spin spaces. The von Neumann entropy is expressed 
in terms of pa as 

5 = -trp^lnp^. (42) 

Therefore, we shall henceforth use the matrix Z instead of M. Our strategy is to find 
a relation between Z matrices for different size systems. Once we find the relation, we 
can construct Z matrices recursively and obtain the von Neumann entropy exactly. 

5.1. VBS states on square ladders 

5.1.1. 2-leg square ladder Let us first study the VBS state on the ladder shown in Fig. 
[9l In this case, the overlap matrix is proportional to 
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(43) 



where n = denotes the horizontal length of A and boundary sites 1 and 2 are shown 
in Fig. 131 The product in the above equation can be expanded as 

n ii-n^-n,)= n (44) 

where the sum runs over all subsets F of B^. Since the measure of integration is invariant 
under the local change of variable, Cli — t- —Cli, we see that the integrand in Eq. ( l43|) must 
contain even number of (ij's. (We give a graphical and combinatorial explanation of this 



fact in terms of self-avoiding loops and strands in Appendix A ) Therefore, Z„(si,S2) 
must be the following form: 

Zn{si, S2) = a„ - b„si ■ S2. (45) 

Then the 4x4 matrix Z„((Ti, 0^2) is written in terms a„ and b„ as 

/ a„ - b„ \ 

a„ + b„ -2b„ 

-2b„ a„ + b„ 

a„ - b„ y 

where a„ and b„ are coefficients depending on n. The eigenvalues of the matrix Z„((Ti, 0^2) 
are a„ + 3b„ and 3-fold degenerate a„ — b„. We shall now compute Z„ by induction on 
n. The relation between Z^+i and is given by 



v 



(46) 



Z„,+i(ai, (T2) 





dfl2 


Air 


Air 


dtli 


dVL2 


Air 


Air 


dtli 


dCl2 


Air 


Air 



l + Cli- ai)(l + ^2 ■ a2)(l - (li ■ ^2) Zni^i, (12) 
l + Cli- ai)(l + (12 ■ (72) (1 - ^1 ■ f^2)(a„ - bndi ■ ^2 



a„ + b„,)(fii ■ CTi)(fi2 ■ 0^2) (^1 ■ ^2) 



(47) 

Here we have used the fact that the integration over odd number of Oj's is zero (see 



Appendix A[ ). To perform the integration over Qi and ^2; the following relation shown 

(48) 



in the lemma 3.3 in Ref. [IE] is useful: 

da 



Air 



(si ■ fii)(fii ■ S2) = qsi ■ S2, 



where g = 1/3. Using this formula, we obtain 

Z„+i(ai, ^2) = a„,+i - K+i^i ■ (T2 = (a„, + gb„) - g^(a„ + b„)ai ■ (T2- (49) 
Therefore, the recursion relation between {a„+i, b^i+i} and {a„, b„} is given by 

3n+l \ /I q \ I Sr 



■>n+l 



q^ 



(50) 
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Here, the initial values ao and bo are 1 and 0, respectively. Equation (1501) can be easily 
solved and a„ and b„ for arbitrary n are given by 

an = ^[4« - + v^(^: + ^")], b„ = -^(^^ - z!!), (51) 



with ^-1- = (5 ± vl9)/9. According to Eq. (14T!) . the eigenvalues of the reduced density 
matrix are given by 



= P3(n) = = (i + 3j]2 + ;(i_,^). ^ (53) 

where x„ = b„/a„. From the above explicit expressions, one can calculate the von 
Neumann entropy for any n. The exact results are totally consistent with those obtained 
from the numerical approach in the previous section (see Table |3]). In the large-block-size 
limit, — )■ 1/(4 + vT9). In this limit, the von Neumann entropy is obtained as 

= -pi(oo) Inpi(oo) - 3p2(oo) lnp2(oo) = 1.2988696. (54) 

Therefore, the von Neumann entropy per valence bond is given by S'^^ /2 = 0.6494348, 
which is strictly less than ln2 = 0.6931471... even in the limit of A^^^ oo. 





N, = l 


N, = 2 


N, = 3 


N, = i 


N, = 5 




= 2 


Exact 
MC 


0.6553433 
0.6553431 


0.6498531 
0.6498533 


0.6494635 
0.6494621 


0.6494368 
0.6494342 


0.6494349 
0.6494373 


Ny 


= 3 


Exact 
MC 


0.6413153 
0.6413145 


0.6325619 
0.6325626 


0.6316999 
0.6316999 


0.6316095 
0.6316080 


0.6315995 
0.6315866 



Table 3. 5/|A^| with \Aa\ = Ny for various obtained by algebraic (Exact) and 
numerical (MC) methods. The numerical values are rounded off to 8 decimal places. 



5.1.2. 3-leg square ladder Let us next study the VBS state on the 3- leg ladder shown 
in Fig. [TOl In this case, the overlap matrix is proportional to 

^n(ai,a2,^3) = / (n^) (l + ^l-'^l)(l + ^2-52)(l + fi3-'?3) J] (l-^^'^i), 

where the boundary sites 1, 2, and 3 are shown in Fig. [TDl Let us set 

Zn{si, 5*2, S3) = a„ — b„si ■ S2 — c„s*2 ■ S3 + d„si ■ S3. (55) 
Note that ao = 1, bo = Cq = do = and 

Zn{^l, <?2, O's) = a„ - b„(Ti ■ (T2 - Cn^2 ' <?3 + dn<?l ' (?3 (56) 
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Figure 10. VBS state on a 3-leg square ladder. 



is 8 X 8 matrix. We shall now compute Z„ by induction on n. The relation between 
Zn+i and Zn is given by 

3 



,i=l 



I (l+fii-<Ti)(l + fi2-a2)(l + 1^3-^3) 



X (1 - fii ■ Q2){1 - ^2 ■ ^3)Zn{ni, ^2, ^s)- (57) 

Now we will use Eq. fH5]) . At the beginning, we exclude from the formula above the 
vector Cli, then ^3, and finally As a result, we will have the recursion relation 
between {a„+i, b„+i, c„+i, d„,+i} and {a„, b„, c„, d„} as 



/ 1 g q \ 



g2 g,2 g,3 g3 
g2 g,2 g,3 g3 
y g3 q3 q,3 q,2 y 



/ an. \ 



d, 



bn+1 
Cn+1 

where g = 1/3. Recalling that ao = 1, bo = Cq 
hn = Cn for all n > 0. Therefore, one obtains 

/ a„ , 1 \ / 1 2q g2 \ / a„ \ 



(5J 



0, we immediately note that 



3n+l 
bn+1 

d,. 



g^ 



g^ g3 



(59) 

y g" 2g'' 

It is now straightforward to obtain the coefficients a^, b^, and d^ by using the above 
recursion relation. 

We shall now calculate the eigenvalues of the reduced density matrix and obtain 
the von Neumann entropies. One can confirm that the eigenvalues of the 8x8 matrix 
Zn in Eq. fl56l) are given by 

Al,n = ^2,n = — 3d„, (60) 



^3,n — A4 „ — A5 „ — Ag 



2bn + d^ 



A7,n — A 



8,n 



4b. 



dn- 



(61) 
(62) 
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Therefore, the eigenvalues of the reduced density matrix are given by 

Mn) = p,in) = 2(1 _ 3^„)2 + 4(1 - 2x^ + yl)^ + 2(1 + 4x„ + y^)^ ' ^^^^ 

/ \ _ _ ( ) — (1 — 2Xn + Vn)^ 

= = 2(1 - 3y,y + 4(1 - 2(1 +4.„ + y„r <«^) 

where 

Xn = — , yn = — • (66) 
Then the von Neumann entropy is expressed as 

8 

Sf,''^ = -J2Mn)\np,{n). (67) 

1=1 

The results are summarized in Table [31 One can again see that the results obtained by 
Monte Carlo integration are consistent with the exact results. Finally, let us consider 
the large-block-size limit. Since the matrix appearing in the recursion relation f l59|l is a 
positive and irreducible matrix, the largest eigenvalue is unique. Let vpp = {vi,V2,V3)'^ 
be the corresponding Perron- Frobenius vector. Note that denotes matrix transpose. 
Then Xn and ?/„, in the limit of n — >■ oo are given by 

X = lim Xn = — , y = hm ?/„ = — . (68) 

n—^oo y-^ n— >oo 

From the numerically obtained vpp, numerical values for x and y are given by 

X = 0.1203998879..., y = 0.0471631199.... (69) 
and hence the von Neumann entropy in the large-block-size limit is obtained as 

8 

S^''^ = - ^pi(oo) Inpi(oo) = 1.8947948. (70) 

i=l 

Therefore, the von Neumann entropy per valence bond is = 0.6315983, which is 

again less than In 2. 

5.2. VBS states on hexagonal ladders 

5.2.1. 2-leg hexagonal ladder Let us now study the VBS state on the hexagonal lattice 
shown in Fig. [TTl which corresponds to Ny = 4. 

Along the same lines as the square-ladder case, the Z matrix is written as 

Zni^l, ^2) = an + bn(Ti ■ 0^2, (71) 

and the recursion relation between {a„+i, b„+i} and {a„, b„} is given by 



an+i \ / 1 

bn+i / \ q'^ j \ K 



(72) 
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Figure 11. VBS state on a 2-leg hexagonal ladder. It is cut by the reflection line 
indicated by the broken line. 



where g = 1/3 and the initial values ao and bo are 1 and 0, respectively. The above 
recursion relation can be solved and the coefficients a„ and b„ for arbitrary n is given 
by 

^ ,mzl-zl) + VlQ^{zl + z-)l b„ = -^(2:-^!!)(73) 



2y/mn ^ ^ 2v^T627 



with z± = (41 ± a/1627)/81. Since the eigenvalues of Zn in Eq. ( I7T1) are a„ — 3b„, and 
3- fold degenerate a„ + b„, the eigenvalues of the reduced density matrix is given by 

= (i-3i.!)^+3aV.„)- ^ ('^' 

p.(„) = p,{n) = p4n) = (75) 

where x„ = b„/a„. The results of the von Neumann entropy are summarized in Table 
HJ In the large-block- size limit, x„ approaches to the value 

X = lim Xn = , (76) 

and hence we obtain the von Neumann entropy 

4 

^^'^s = - J]pi(oo) \npi{oo) = 1.3781854. (77) 

1=1 

The von Neumann entropy per valence bond is given by = 0.6890927, which is 

strictly less than In 2. 



5.2.2. 3-leg and hexagonal ladder Next we shall study the VBS states on the 

3-leg and 4-leg hexagonal ladders. As an illustration, in Fig. [121 we show an example 
of 3-leg hexagonal ladders. We ffist study the 3-leg ladders. In this case, the Z matrix 
is defined similarly to Eq. ( !56l) : 

^n(<?l, 0^2, (Ts) = a„ + b„(Ti ■ (T2 + C„a2 " (?3 + d„(Ji ■ (T3, (78) 
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N, = l 


N, = 2 


N, = 3 


N, = i 


N, = 5 


Ny=A 


Exact 
MC 


0.6891577 
0.6891575 


0.6890932 
0.6890924 


0.6890927 
0.6890929 


0.6890927 
0.6890925 


0.6890927 
0.6890840 


Ny = Q 


Exact 
MC 


0.6878024 
0.6878027 


0.6876554 
0.6876558 


0.6876523 
0.6876513 


0.6876522 
0.6876537 


0.6876522 
0.6875899 


Ny = 8 


Exact 
MC 


0.6871245 
0.6871243 


0.6869344 
0.6869385 


0.6869295 
0.6869363 


0.6869293 
0.6868834 


0.6869293 
0.6867750 



Table 4. 5/|Ayi| with |A^| = Ny/2 for various Nx obtained by algebraic (Exact) and 
numerical (MC) methods. The numerical values are rounded off to 8 decimal places. 




with ao = 1 and bo = Cq = do = 0. The recursion relation becomes slightly complicated 
compared to the square cases due to the wiggly nature of the hexagonal ladders. It is 
composed of 2 steps as follows: 



(32m+l) b2m+l; ^2m+l, d2m+l)'^ 
(32m+25 b2m+25 C2m+25 d2m+2)"'" 

where two matrices Ti and T2 are given by 



TiT2(a2m-l; ^2m-l, '^2m-l, d2m-l) 
T2ri(a 2m, ^2m, ^2m, d2m) , 



(79) 
(80) 



/ 1 



q^ 



\ 



q q q q 



6 



q q 



4 4 

q q 

6 4 

q q 



( 1 



To 



q q q 



q q q q 

g3 q,5 g4 g,6 



(81) 



y 

respectively, with q = 1/3. Using the above recursion relation, one can obtain 
the coefficients {a„, b„, Cn, d.„}. By diagonalizing the matrix in Eq. (178|) with these 
coefficients, we can calculate the von Neumann entropy along the same lines as Sec. 
15.1.21 The obtained results are summarized in Table HI which are compared with 
numerical results obtained in Sec. 14.21 We now consider the large-block-size limit. 
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Let vpp be the Perron- Frobenius vector of T1T2. This vector is numerically given by 

vpF = (1, 0.03734899, 0.03734899, 0.0046503105)'^ (82) 
Using the above values, we obtain the von Neumann entropy in the large-block-size limit 



as 



5£i^g = 2.0629566, 



(83) 



and hence the von Neumann entropy per valence bond is given by = 0.6876522, 

which is strictly less than In 2. 

As is obvious, the argument so far can be straightforwardly generalized to m-leg 
ladders with m > 4 |104j . In Table IU we show the results for the 4-leg hexagonal 
ladders. All the analytically obtained results indicate that the boundary correlation 
does not vanish even in the large-block- size limit and support the conjecture in Sec. 
[3l i.e., the eigenvalues of the reduced density matrix can be different. A numerical 
implementation of the transfer matrix technique for the m-leg ladders with m > 4 will 
be discussed elsewhere |104] . 



6. Vertical ladders 

So far, we have studied square or hexagonal ladders in a horizontal geometry. In this 
section, we shall study the opposite limit, i.e., the ladders in a vertical geometry. Let 
us first consider the square ladder with A^^, = 1 and the height Ny shown in Fig. [TST a). 



(a) 



A 

4 



>^ Ny 



T 
T 



B 

€) 

T 

T 

T 

T 
€) 




Figure 13. (a) Square lattice vertical ladder, (b) Hexagonal lattice vertical ladder. 



Then the Z matrix is given by 



Z{Ny]ai,...,aNy) 



\i=i J j=i k=i 



(84) 
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Similarly to the horizontal ladder cases, we can use the integral formula Eq. f l48p and 
rewrite Eq. fl8^ as 

Z{Ny) = J2 q'^i-qr'-'K-qr'-'' ■ ■ ■ (-g)^--'="-n^n-^^J(^.3-c?M) • • • (^..„-.-^..J, (85) 

{«!,..., j2n} 

where 1 < ii < 12 < ■ ■ ■ < i2n-i < Hn < Ny. Here we have abbreviated 
Z{Ny) = Z{Ny;ai, ...,aNy)- The above formula can also be obtained using loops and 
strands shown in Appendix A Note that Z{1) is 2 x 2 identity matrix. We now look 
for a recursion relation between Z{Ny + 1) and Z{Ny). We first decompose Z{Ny) into 
four matrices as 

3 

Z{Ny) = Z''{Ny) + J2z''iNy), (86) 

a=l 

where Z^{Ny) trivially acts on the A^^^-th vector space Vn^ while Z°'{Ny) acts on it as a". 
This decomposition is graphically shown in Fig. [HI Using the graphical representation, 



ZiNy) Z\Ny) 



Z"(Ny) 



Ny 



Ny-l 



3 

«=1 



id 



Figure 14. Decomposition of the matrix N{Ny). 



we find the following set of recursion relations: 

Z%Ny + 1) = Z{Ny) ® id, 

Z^{Ny + 1) = - q^Z{Ny - 1) ® a" ® a" - g 



(87) 



Z'^iNy). I (g)id) 



a 



Z{Ny + l) =Z'{Ny + l) + J2N''iNy + l), (89) 

a=l 

where '.' denotes matrix multiplication while '®' denotes tensor product, and 'id' is 
the 2x2 identity matrix. Solving the above recursions by mathematica, we obtain the 
overlap matrix Z{Ny) exactly up to Ny = 12. Then we numerically diagonalize Z{Ny) 
and calculate the von Neumann entropy as a function of Ny. Figure [TH a) shows the 
result. The obtained values and also fitting parameters in Eq. ( 138|) show good agreement 
with the Monte Carlo results. 
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Figure 15. (a) The von Neumann entropy per valence bond as a function of Ny for 
the square lattice VBS state with = 1. The blue curve shows the fitting curve (Eq. 
(pl) ) with C = 0.0819(3), A = 0.91(1), a = 0.6113(3). (b) The von Neumann entropy 
per valence bond as a function of Ny for the hexagonal lattice VBS state with = 1. 
The blue curve shows the fitting curve (Eq. ([38])) with C = 0.008081(5), A = 0.985(1), 
a = 0.685068(5). 



Let us next consider the hexagonal ladder with A'a; = 1 and the height Ny shown 
in Fig. [TST b). It is convenient to set A'j^ = 2m — 1 since Ny is always odd in this case. 
In this case, from Eq. ( HOl) . the Z matrix is given by 



Z{m]ai, ...,am) 



■N, 



(90) 

yi=i " / j=i k=i 

Along the same lines as the square ladder case, we obtain the set of recursion relations 
for the Z matrix: 

Z°(m + 1) 

^ m— 1 
(g)id 



Z(m) (g) id, 
Z"(m + 1) = q^Z{m - 1) ® a" ® a" + 



Z'^im) 



a 



(91) 
(92) 



Z{m + 1) = Z°(m + 1) + ^ Z"(m + 1). 



(93) 



a=l 



Here we have abbreviated Z{m) = Z{m\ (Xi, (?,„). Again, solving the above recursions, 
we find the Z matrix up to m = 12 {Ny = 23). The numerically calculated von Neumann 
entropy by this method is shown in Fig. [T^ b). The obtained values and the fitting 
parameters also show good agreement with the Monte Carlo results. The obtained 
results for square and hexagonal ladders with N^ = 1 strongly support that the assumed 
fitting function Eq. (138|) correctly describes the behavior of large- A'^y systems. 
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7. Conclusion 

In conclusion, we have studied entanglement properties of the VBS states on two- 
dimensional graphs with reflection symmetry. We have shown that the reflection 
symmetry permits us to develop an eflicient method to study the reduced density 
matrix of the subsystem which is a mirror image of the other one. We found the 
relation between the reduced density matrix and the overlap matrix whose element is 
given by the inner product between degenerate ground states of the block Hamiltonian 
defined in the subsystem. This relation enables us to conjecture that the von Neumann 
entropy per valence bond on the boundary between the subsystems is strictly less than 
In 2 in 2d systems. This conjecture means that the entanglement spectrum is not flat 
in this system and is in contrast to the case of Id VBS where the reduced density 
matrix is proportional to the identity matrix in the large-block limit. We have given 
a holographic interpretation of the reduced density matrix in terms of the spin chain 
along the boundary between the subsystems, which describes a hidden correlation among 
the degenerate ground states of the block Hamiltonian. To conflrm the conjecture, we 
have numerically studied the eigenvalues of the reduced density matrix for flnite square 
and hexagonal graphs by combining Monte Carlo integration with exact diagonalization. 
From the eigenvalues of the reduced density matrix, we have calculated the von Neumann 
entropy and have found that the von Neumann entropy per valence bond is well fltted 
by a constant term with a function decaying algebraically in the length of the boundary. 
The constant is strictly less than In 2 and supports our conjecture. We have also 
analytically and algebraically studied the quasi-ld cases where graphs are on ladders. 
The analytical results are totally consistent with the numerical results and strongly 
support our conjecture. Although our conjecture is very plausible, we still lack a rigorous 
proof of it in the inflnite 2d system. Therefore, it would of course be interesting to 
study the problem by a completely different approach based on inequalities and perhaps 
reflection positivity. 
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Appendix A. Graphical interpretation of Z matrix 

In this appendix, we iUustrate a graphical interpretation of the overlap matrix Z (see 
Eq. f H3|l ) in terms of loops and strands. As we have already mentioned in Sec. 15. ![ the 
terms containing odd number of f2j's in 

J] (i-a-fi,)=5^ (A.i) 

vanish after integration over f2j since the measure of integration is invariant under the 
local change of variable, fij — )■ — fij. This can be graphically explained as follows. 
Suppose that one vertex z G A is connected to two vertices. They are labeled 1 and 2 
(see Fig. lAl( a)). Then we consider the following integral: 

I _ ti., ■ fii)(i -^^■^2) = I ^[1 + (a ■ n,){n, ■ n^)]- (a.2) 

From the lemma 3.3 in Ref. [H], we can show 

dn,{n, ■ ■ n^) = —n^ ■ ^2 (a.3) 

and hence we obtain 

//-70 
^(1 - a ■ ^^i)(l -Vti-Q^) = l + q^i- ^2, (A.4) 

where g = 1/3. This procedure is graphically shown in Fig. lAl( a). Next, suppose that 
one vertex i is surrounded by three vertices labeled by 1, 2, and 3. Then we encounter 
with the following integral: 

j ^(1 - a ■ - ■ n2){i - ■ n^)- (a.s) 

Using Eq. ( 1A.3I) . again, we obtain 

^(1 - ■ n,){i - a, ■ n2){i - a ■ ^3) 

= 1 + ■ ^2 + g^2 -^3 + g^3 ■ ^1- (A.6) 

This procedure can also be graphically interpreted as shown in Fig. lAl( b). Therefore, 
as we have seen, the integrand containing odd number of fij's will vanish. In terms of 
the graphical representation, the matrix Z„ for the 2-leg square ladder can be written 

as 

Z4a,,a2) = ^g^B{C)^-A^.{C) _j2q^-^NMc'})q-NMc'})^^^ . ^^)^ (a.7) 

{C} {C} 

where A''b(C) and Ni^{C) denote the numbers of colored edges (bonds) and loops 
in the configuration C, respectively. The minus sign stems from the product like 
(1 + f2i ■ ai)(l — Qi ■ n2)(l + ^2 ■ ^2)- In the first sum over {C}, only the closed 
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1/2 1 21 2 



(b) 




Figure Al. (a) Graphical representation for Eq. (jA.4[) . The mark (g) denotes the 
integration over the spherical angle fli. The broken and colored (thick) lines have 
weight 1 and q = 1/3, respectively, (b) Graphical representation for Eq. (jA.6[) . 
Notations are the same as (a). 



n = l 



n = 2 



n = 3 



□ 



c 



i:t:i 



Figure A2. Closed loop C and open strand C configurations for the length n = 1,2, 
and 3. The weight of each configuration is assigned at the bottom right of each figure. 



loops are allowed. On the other hand, in the second sum over {€'}, open ends from 1 
and 2 are allowed (see Fig. IA2I) . Therefore, we can obtain the coefficients a„ and b„ in 
Eq. fH5l) by a completely graphical or combinatorial way. From Fig. IA2t one can read 
ai = 1, 32 = 1 + g^ as = 1 + 2q^ + and bi = ba = + g^, = + + + q^. 
Although we only show the case of 2-leg square ladders for simplicity, one can easily 
generalize this approach to any AKLT ladders. 
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